% This m-scripts checks the rank of the coefficient matrix compiled by
% stacking structural restrictions across mutliple environments indexed by
% 's' in LQT(2019), e.g. classes with different sizes. 

% Exclusion restriction: individual and contextual effects do not vary
% across environments, but peer effects do. 

% This m-scripts shows that the rank of such a coefficient matrix
% is generically one fewer than # of unknown structural parameters. 

% Date: Nov 13, 2018

clear; clc; 

%% define parameter values
K = 4;  % # of regressors
beta = rand(K,1); gamma = rand(K,1); gamma(1) = 0; % individual and contextual effects
% peer effect for small and large classes
lam = [.7; .8]; % rand(2,1); % [.65; .75]; % rand(L,1); % [.7 ; .8; .6];   
L = length(lam); % # of environments

%% solve for (a_k,b_k) for k < K in each environment indexed by l<K
abMat = zeros(2,K-1,L);   
% Dimensions of 'abMat'
% - 1st dim: (a,b) 
% - 2nd dim: index for regressors (k) 
% - 3rd dim: class size, small or large
M = zeros(K,L);
% Dimensions of 'M'
% - 1st dim: index for regressors (k) 
% - 2nd dim: class size, small or large
for i = 1:L, % loop in environments (e.g. class size)
    for k = 1:K-1,
        abMat(:,k,i) = [beta(k) beta(K); gamma(k) gamma(K)]\[1;-lam(i)];
        M(k,i) = (beta(k)+gamma(k))/(1-lam(i));
    end;
    M(K,i) = (beta(K)+gamma(K))/(1-lam(i));
end;
a = squeeze(abMat(1,:,:)); % 'a' and 'b' are both (K-1)-by-L matrices
b = squeeze(abMat(2,:,:)); 
%% construct the coefficient matrix and the vector of constants for the
% linear system after stacking L environments
Pi = zeros( L*(2*(K-1)+K) , L+2*K ); 
Mat1 = []; Mat2 = []; D = [];
for ci =  1:L, % loop in environments
    % fill one out of the first L columns in Phi
    vec = zeros(L,1); vec(ci) = 1; 
    vci = zeros( 2*(K-1)+K,1 ); 
    vci( 2*(K-1)+1:end ) = M(:,ci); 
    vci(K:2*(K-1)) = ones(K-1,1);
    Pi(:,ci) = kron(vec,vci);
    % construct component matrices for the other 2K columns
    mat1 = zeros(2*(K-1)+K,K); 
    mat1(1:K-1,:) = [diag(a(:,ci)) b(:,ci)];
    mat1( 2*(K-1)+1:end ,:) = eye(K);  
    Mat1 = [Mat1;mat1];     
    mat2 = zeros(2*(K-1)+K,K); 
    mat2(K:2*(K-1),:) = [diag(a(:,ci)) b(:,ci)];    
    mat2( 2*(K-1)+1:end ,:) = eye(K);    
    Mat2 = [Mat2;mat2];    
    cvec = [ones(K-1,1); zeros(K-1,1); M(:,ci)];
    D   = [D; cvec];
end;
Pi(:,L+1:end) = [Mat1 Mat2];

%% report the ranks and check for bugs in coding up the system
disp(sprintf('Coef Matrix using L = %g environments with K = %g:',L,K));
disp('(# rows, # columns, # rank)');
disp([size(Pi), rank(Pi)]); 
max(abs(Pi*[lam;beta;gamma] - D)); % check for bugs in linear system

%% add a single exclusion restriction and check ranks
disp('*************************************');
disp('Add a linear restriction: 1st regressor has no contextual effect');
R = [zeros(1,L+K),1,zeros(1,K-1)]; 
max(abs([Pi;R]*[lam;beta;gamma] - [D;0])); % check for bugs in coding
disp('(# rows, # columns, # rank)');
disp([size([Pi;R]), rank([Pi;R])]); 

